* Analysis of Double List Experiment (Reproduces Figure 1)
						
* Opens data
use "$directory/List experiment/data/CRON2W6_PT_e01_translation.dta", clear


** Makes a treatment indicator for List A:

gen D1 = .
replace D1 = 1 if w6ptq58b < 6						// Treatment group ("the government should be able to ignore court rulings that are regarded as politically biased." )
replace D1 = 0 if w6ptq58a < 6						// Control group -- not shown sensitive statement

* Item-count outcome #1: 
recode w6ptq58a (6 9 = .), gen(w6ptq58a_2)
recode w6ptq58b (6 9 = .), gen(w6ptq58b_2)

gen count1 =. 
replace count1 = w6ptq58a_2 if D1 == 0
replace count1 = w6ptq58b_2 if D1 == 1


** Makes a treatment indicator for List B:
gen D2 = .
replace D2 = 1 if w6ptq59a < 6					// Treatment group ("the government should be able to ignore court rulings that are regarded as politically biased." )
replace D2 = 0 if w6ptq59b < 6					// Control group -- not shown sensitive statement


* Item-count outcome #2: 
recode w6ptq59b (6 9 = .), gen(w6ptq59b_2)
recode w6ptq59a (6 9 = .), gen(w6ptq59a_2)

gen count2 =. 
replace count2 = w6ptq59a_2 if D2 == 1
replace count2 = w6ptq59b_2 if D2 == 0


** Direct question: 
recode JUDCNSTR2 (9 =.) (1 2 = 1 "agree") ///
	(3 4 5 = 0 "disagree"), gen(direct)

** Double-list analysis using the "KICT" package (see Tsai 2019)
gen T = . 										// Create a dummy, T, where = 0 those in short-list first and long-list second and 1 = opposite
replace T = 1 if D1 == 1 & D2 == 0				// respondents who received the short-list in the second item count
replace T = 0 if D1 == 0 & D2 == 1				// respondents who received the short-list in the first item count	
	
	
** Comparing list estimate to direct estimate using equality of proportions tests (results presented in footnote 11):
sum direct 												// n = 365, mean = .21 
kict ls count1 count2, ///
condition(T) nnonkey(4 4) estimator(linear) duallist  	// n = 361, mean = .177

prtesti 361 .177 365 .210, level(95)	


* ------------------ *
* 		FIGURE 1	 *
* ------------------ *

	
* Comparing list estimate to direct using 1000 bootstraps 
set seed 1111

tempfile `compdata'
postfile results delta raw using "$directory/List experiment/data/bootstrap_dif", replace
forvalues i=1/1000{
preserve
bsample

kict ls count1 count2, condition(T) nnonkey(4 4) estimator(linear) duallist 
mat A = e(b)
local g=A[1,1]

sum direct 
local r= r(mean)



post results (`g') (`r') 
	
restore
}
postclose results


* Plot results (Fig. 1)
use "$directory/List experiment/data/bootstrap_dif.dta", clear


* Makes a variable for difference in estimates across bootstrap replications for direct and dle estimates: 
gen dif = delta - raw		
su dif  					// mean = -.0336376, sd = .0606777. SD in this bootstrapped sampling distribution is estimate for standard error of difference in prevalence estimates.  
hist dif, xline(0)
count if dif > 0			// number of estimated differences in prevalence that fall above zero is estimated p-value for difference in prevalence


mean dif 
ereturn list
mat A = e(sd)*1.96
mat B = e(b)+A
mat C = e(b)-A
mat D = e(b)

mean delta
ereturn list
mat E = e(sd)*1.96
mat F = e(b)+E
mat G = e(b)-E
mat H = e(b)

mean raw
ereturn list
mat I = e(sd)*1.96
mat J = e(b)+I
mat K = e(b)-I
mat L = e(b)


mat P = 1,.7,H[1,1],G[1,1],F[1,1]\ 	///
		2,1,L[1,1],K[1,1],J[1,1]\ 	///
		3,.4,D[1,1],C[1,1],B[1,1]
		
		
* Plot estimates and differences in estimates
frame create plot1
cwf plot1
svmat P
rename P1 outcome
rename P2 Ypos
rename P3 b
rename P4 ll
rename P5 ul

tw ///
	(rspike ll ul Ypos if outcome == 1, horizontal lcolor(black)) ///
	(rspike ll ul Ypos if outcome == 2, horizontal lcolor(black)) ///
	(rspike ll ul Ypos if outcome == 3, horizontal lcolor("255 64 64")) ///
	(scatter Ypos b if outcome == 1, msym(O) mc(black) mfcolor(white) msize(medlarge) mlab(b) mlabformat(%5.2fc) mlabposition(12) mlabcolor(black) mlabsize(med)) ///
	(scatter Ypos b if outcome == 2, msym(O) mc(black) mfcolor(white) msize(medlarge) mlab(b) mlabformat(%5.2fc) mlabposition(12) mlabcolor(black) mlabsize(med)) ///
	(scatter Ypos b if outcome == 3, msym(O) mc("255 64 64") msize(medlarge) mlab(b) mlabformat(%5.2fc) mlabposition(12) mlabcolor("255 64 64") mlabsize(med)), ///
	yscale(range(.1 1.2) lc(none)) ytit("") ///
	ylab(.7 `" "List experiment" "prevalence estimate" "' 1 `" "Direct question" "prevalence estimate" "' .4 `" "Difference in" "estimated prevalence" "', nogrid) ///
	xtit("") xlab(-.3 -.2 -.1 0 .1 .2 .3, notick nogrid) xline(0) xscale(lc(none)) ///
	legend(off) plotregion(lcolor(black) lwidth(medium)) ///
	tit("") ///
	name(Figure1, replace)

graph display, xsize(4.5) ysize(2) scale(1.4)
graph export "$directory/Output/Figure1.pdf", as(pdf) replace

